Simulation method

ABSTRACT

(a) Regarding a particle system S in which the number of particles is N, the mass of each particle is m, and inter-particle interaction potential energy can be expressed by εf, α greater than 1, γ equal to or greater than 0 and equal to or smaller than d, and δ equal to or greater than 0 are determined using a dimension number d of a space where the particle system S is arranged to obtain the number N′ of renormalized particles by N′=N/α d , to obtain the mass m′ of each of the renormalized particles by m′=mα δ /α γ , and to obtain a renormalized interaction coefficient ε′ by ε′=εα γ . (b) Molecular dynamics calculation is carried out on a particle system S′ in which the number of renormalized particles is N′, the mass of each particle is the mass m′ of each renormalized particle, and inter-particle interaction potential energy is expressed by ε′f.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation application of PCT Patent ApplicationPCT/JP2009/005747 filed on Oct. 29, 2009, which is based on and claimspriority of Japanese Patent Application JP2008-324090 filed on Dec. 19,2008, the entire contents of which are incorporated herein by reference.

BACKGROUND OF THE INVENTION

A) Field of the Invention

The present invention relates to a simulation method and a program andin particular, to a simulation method using molecular dynamics and aprogram which causes a computer to execute the simulation method.

B) Description of the Related Art

Computer simulations using molecular dynamics are being carried out. Inthe molecular dynamics, the equation of motion of particles constitutinga system serving as a simulation target is analyzed numerically. If thenumber of particles constituting a system serving as a simulation targetincreases, the amount of necessary calculation increases. With thearithmetic capacity of an existing computer, it is only possible tocarry out the simulation of a system having particles of about a hundredthousand.

In order to reduce the amount of calculation necessary for thesimulation, for example, as described in NPL 1 (NPL 1: Inamura,Takezawa, and Shamoto, “On Variable Scale Molecular Dynamics SimulationBased on Renormalization Technique”, Transactions of the Japan Societyof Mechanical Engineers (A), 1997, Vol. 63, No. 608, p. 202-207),attempts have been made to apply a renormalization technique to themolecular dynamics.

The invention relating to this application has been made by theinventors and described in PTL 1 (PTL 1: JP-A-2006-285866).

SUMMARY OF THE INVENTION

As described in “3. Numerical Calculation Method” of NPL 1, in thesimulation method of NPL 1, it may be impossible to evaluate atemperature explicitly.

One object of the invention is to provide a novel simulation methodusing molecular dynamics and a program which causes a computer toexecute the method.

Another object of the invention is to provide a simulation methodcapable of evaluating a temperature explicitly on the basis of arenormalization technique and a program which causes a computer toexecute the method.

An aspect of the invention provides a simulation method. The simulationmethod includes the steps of (a) with regard to a particle system S inwhich the number of particles is N, the mass of each particle is m, andinter-particle interaction potential energy can be expressed by aproduct εf of a dimensionless function f representing an inter-particledistance dependence and an interaction coefficient ε, determining afirst renormalization factor α greater than 1, a second renormalizationfactor γ equal to or greater than 0 and equal to or smaller than d, anda third renormalization factor δ equal to or greater than 0 using adimension number d of a space where the particle system S is arranged toobtain the number N′ of renormalized particles by a conversion equationN′=N/α^(d), to obtain the mass m′ of each renormalized particle by aconversion equation m′=mα^(δ)/α^(γ), and to obtain a renormalizedinteraction coefficient ε′ by a conversion equation ε′=εα^(γ), and (b)carrying out molecular dynamics calculation on a particle system S′ inwhich the number of renormalized particles is N′, the mass of eachparticle is the mass m′ of each renormalized particle, andinter-particle interaction potential energy is expressed by a productε′f of the dimensionless function f and the renormalized interactioncoefficient ε′.

Another aspect of the invention provides a simulation method. Thesimulation method includes the steps of (a) with regard to a particlesystem S in which the number of particles is N, the mass of eachparticle is m, and inter-particle interaction potential energy can beexpressed by a product εf of a dimensionless function f representing aninter-particle distance dependence and an interaction coefficient ε,determining a first renormalization factor α greater than 1 and a thirdrenormalization factor δ equal to or greater than 0 using a dimensionnumber d of a space where the particle system S is arranged to obtainthe number N′ of renormalized particles by conversion equationN′=N/α^(d) and to obtain the mass m′ of each renormalized particle by aconversion equation m′=mα^(δ), and (b) carrying out molecular dynamicscalculation on a particle system S′ in which the number of renormalizedparticles is N′, the mass of each particle is the mass m′ of eachrenormalized particle, and inter-particle interaction potential energyis expressed by the product εf of the dimensionless function f and theinteraction coefficient ε.

Yet another aspect of the invention provides a simulation method. Thesimulation method includes the steps of (d) with regard to a particlesystem S in which the number of particles is N, the mass of eachparticle is m, and inter-particle interaction potential energy can beexpressed by a product εf of a dimensionless function f representinginter-particle distance dependence and an interaction coefficient ε, andan XYZ orthogonal coordinate system in a space where the particle systemS is arranged, on the basis of a renormalization factor α_(X) in an Xdirection, a renormalization factor α_(Y) in a Y direction, and arenormalization factor α_(Z) in a Z direction greater than 1, obtainingthe number N′ of renormalized particles by a conversion equationN′=N/(α_(X)α_(Y)α_(Z)), and further determining a renormalization factorδ equal to or greater than 0 to obtain the mass m_(X)′ of eachrenormalized particle regarding the X direction by a conversion equationm_(X)′=mα_(X) ^(δ), to obtain the mass m_(Y)′ of each renormalizedparticle regarding the Y direction by a conversion equationm_(Y)′=mα_(Y) ^(δ), and to obtain the mass m_(Z)′ of each renormalizedparticle regarding the Z direction by a conversion equationm_(Z)′=mα_(Z) ^(δ), and (e) carrying out molecular dynamics calculationon a particle system S′, in which the number of renormalized particlesis N′, the mass of each particle is the mass m′ of each renormalizedparticle, and inter-particle interaction potential energy is expressedby the product εf of the dimensionless function f and the interactioncoefficient ε, using the mass m_(X)′ of each renormalized particleregarding the X direction for an equation of motion in the X direction,the mass m_(Y)′ of each renormalized particle regarding the Y directionfor an equation of motion in the Y direction, and the mass m_(Z)′ ofeach renormalized particle regarding the Z direction for an equation ofmotion in the Z direction.

The number N′ of renormalized particles is smaller than the number N ofparticles. For this reason, molecular dynamics calculation can becarried out on the particle system S′ with a small amount of calculationcompared to a case where molecular dynamics calculation is carried outon the particle system S.

The renormalization factor δ equal to or greater than 0 is arbitrarilydetermined and calculation is carried out, thus it is possible to carryout a simulation flexibly.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a diagram illustrating a method of deriving a scaleconversion rule.

FIG. 1B is a diagram illustrating a method of deriving a scaleconversion rule.

FIG. 2 is a graph showing a simulation result on aluminum usingrenormalization group molecular dynamics.

FIG. 3A is a diagram illustrating an S-W potential.

FIG. 3B is a graph showing a simulation result on silicon usingrenormalization group molecular dynamics.

FIG. 4 is a graph showing a simulation result when a dispersionrelationship is obtained using renormalization group molecular dynamics.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

First, molecular dynamics (MD) will be simply described. A particlesystem in which the number of particles (for example, atoms) is N, andthe Hamiltonian H is expressed by the following equation is taken intoconsideration.

$\begin{matrix}{H = {\sum\limits_{j = 1}^{N}\left\lbrack {\frac{p_{j}^{2}}{2m} + {\sum\limits_{i = {j + 1}}^{N}{\phi\left( {q_{i} - q_{j}} \right)}}} \right\rbrack}} & (1)\end{matrix}$

Here, the mass of each particle is m, the momentum vector of a particlej is p_(j), the position vector (position coordinate) of the particle jis q_(j), and inter-particle interaction potential energy is φ.

If Hamiltonian H is substituted in a Hamiltonian canonical equation, thefollowing equations of motion for each particle are obtained.

$\begin{matrix}{\frac{\mathbb{d}p_{j}}{\mathbb{d}t} = {- {\sum\limits_{i \neq j}^{N - 1}{\frac{\partial{\phi\left( {q_{i} - q_{j}} \right)}}{\partial q_{j}}\mspace{14mu}\left( {{j = 1},\ldots\mspace{14mu},N} \right)}}}} & (2) \\{\frac{\mathbb{d}q_{j}}{\mathbb{d}t} = {\frac{p_{j}}{m} = {v_{j}\mspace{14mu}\left( {{j = 1},\ldots\mspace{14mu},N} \right)}}} & (3)\end{matrix}$

In the equation (3), v_(j) is the velocity vector of the particle j. Inmolecular dynamics, for each particle constituting a particle system, anequation of motion expressed by the equations (2) and (3) is numericallyintegrated and solved, obtaining the momentum vector (or velocityvector) and the position vector of each particle at each time. With thenumerical integration, in many cases, a Verlet method is used. TheVerlet method is described in, for example, J. M. Thijssen,“Computational Physics”, Cambridge University Press (1999), p. 175.Various physical quantities of a system can be calculated on the basisof the position and speed of each particle obtained by moleculardynamics calculation.

Next, molecular dynamics (this is called renormalization group moleculardynamics) using a renormalization group method will be conceptuallydescribed.

In renormalization group molecular dynamics, a desired system S in whichthe physical quantities will be obtained is associated with a system(this is called a renormalized system S′) having a smaller number ofparticles than the system S. Next, molecular dynamics calculation iscarried out on the renormalized system S′. The calculation result on therenormalized system S′ is associated with the desired system S. Thus, itbecomes possible to calculate the physical quantities of the desiredsystem S with a small amount of calculation compared to a case wheremolecular dynamics calculation is carried out directly on the desiredsystem S. A conversion rule which associates the physical quantity (forexample, number of particles, the mass of each particle, and the like)in the desired system S with the physical quantities in the renormalizedsystem S′ is called a scale conversion rule.

Next, the scale conversion rule which has been found by the inventorswill be described. It is assumed that the desired particle system S inwhich the physical quantities will be calculated has N particles. It isassumed that the entire Hamiltonian H of the particle system S isexpressed as shown in the equation (1).

The Hamiltonian H_(j) of a certain particle j is expressed as follows.

$\begin{matrix}{H_{j} = {\frac{p_{j}^{2}}{2m} + {\sum\limits_{i = {j + 1}}^{N}{\phi\left( {q_{i} - q_{j}} \right)}}}} & (4)\end{matrix}$

For discussion on renormalization, interaction potential energy φ isexpressed as follows by a product of a function f which representsinter-particle distance dependence and is dimensionless and acoefficient ε (this is called an interaction coefficient ε) whichrepresents interaction intensity and has a dimension of energy.

$\begin{matrix}{{\phi\left( {q_{i} - q_{j}} \right)} = {ɛ\;{f\left( \frac{q_{i} - q_{j}}{\sigma} \right)}}} & (5)\end{matrix}$

For example, with regard to an inactive atom, the following interactionpotential is used.

$\begin{matrix}{{\phi(r)} = {ɛ\left\lbrack {\left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6}} \right\rbrack}} & (6)\end{matrix}$

Here, r is an inter-atom distance. The interaction coefficient ε and theconstant σ have values corresponding to the type of atom.

Next, a method of deriving a renormalized Hamiltonian H′ from theHamiltonian H will be described. As described below in detail, therenormalized Hamiltonian H′ is obtained by executing a part ofintegration in a partition function Z(β) for the particle system S,coarsely graining the Hamiltonian, and subsequently rescaling anintegration variable.

For the particle system S, the partition function Z(β) of a canonicalensemble with a constant number of particles is expressed as follows.Z(β)=∫dΓ _(N)exp(−βH(p,q))  (7)

Here, the coefficient β is defined as β=1/(k_(B)T) on the basis of thetemperature T of the system and the Boltzmann constant k_(B). dΓ_(N) isa volume element in a phase space and is specifically expressed asfollows.

$\begin{matrix}{{\mathbb{d}\Gamma_{N}} = {{\frac{1}{W_{N}}{\prod\limits_{j = 1}^{N}\;{{\mathbb{d}p_{j}}{\mathbb{d}q_{j}}}}} = {\frac{1}{W_{N}}D_{p}^{N}D_{q}^{N}}}} & (8)\end{matrix}$

Here, W_(N) is N!h^(3N) (where h is the Planck constant). D_(P) ^(N) andD_(q) ^(N) are respectively defined as follows.

$D_{p}^{N} = {\prod\limits_{j = 1}^{N}\;{\mathbb{d}p_{j}}}$$D_{q}^{N} = {\prod\limits_{j = 1}^{N}\;{dq}_{j}}$

It is difficult to discuss renormalization of a potential other than aharmonic oscillator. Accordingly, it is considered that an interactionpotential having a saddle point is divided into three regions. Since theinter-particle distances r→∞ and r→0 are fixed points, the potential isinvariable at the time of renormalization conversion. Thus, it isassumed that renormalization conversion near the saddle point isestablished for all the inter-particle distances r. Around the saddlepoint is subjected to Taylor expansion, and fluctuations up to asecond-order fluctuation remain. If the position of the saddle point isr₀, and the relative displacement is δu, the first-order term of δu iscleared, and the following equation is obtained.

$\begin{matrix}{{\phi\left( {r_{0} + {\delta\; u}} \right)} = {{\phi\left( r_{0} \right)} + {{\frac{1}{2}\left\lbrack \frac{\partial^{2}{\phi(r)}}{\partial r^{2}} \right\rbrack}_{r = r_{0}}\delta\; u^{2}}}} & (9)\end{matrix}$

Thus, for the particle j, the Hamiltonian which describes near thesaddle point is a follows.

$\begin{matrix}{H_{j} = {\frac{p_{j}^{2}}{2\; m} + {\sum\limits_{i = {j + 1}}^{N}\left\lbrack {{\frac{1}{2}{\phi^{''}\left( r_{0} \right)}{{u_{i} - u_{j}}}^{2}} + {\phi\left( r_{0} \right)}} \right\rbrack}}} & (10)\end{matrix}$

Here, u_(i) and u_(j) respectively represent the displacement from thesaddle points of particles i and j. φ″ is a second order differentialfor the inter-particle distance r.

Next, coarse graining of an inter-particle interaction portion of thepartition function Z(β) will be described. Consideration is first madeon coarse graining in a particle system arranged on a one-dimensionalchain, and is expanded to a particle system arranged on a simple cubiclattice.

As shown in FIG. 1A, at the lattice points on a one-dimensional chain,particles i and j are arranged to be nearest neighbor, and particles jand k are arranged to be nearest neighbor. The particle j is arrangedbetween the particles i and k. The interaction (nearest neighborcoupling) between nearest neighbor particles is indicated by a solidline. It is considered that the particle j is erased, and coarsegraining is carried out. If the interaction in which the particle jinvolves is written, and the displacement u_(j) of the particle jbetween the particles i and k is integrated, for the interactionpotential φ′(q_(i)−q_(k)) of the particles and k after the influence ofthe particle j is renormalized, the following equation is obtained.

$\begin{matrix}\begin{matrix}{{\exp\left( {- {{\beta\phi}^{\prime}\left( {q_{i} - q_{k}} \right)}} \right)} = {\int_{- \infty}^{\infty}\ {{\mathbb{d}u_{j}}{\exp\left\lbrack {{- {{\beta\phi}^{''}\left( r_{0} \right)}}\begin{Bmatrix}{\frac{\left( {u_{i} - u_{j}} \right)^{2}}{2} +} \\{\frac{\left( {u_{j} - u_{k}} \right)^{2}}{2} +} \\{2\frac{\phi\left( r_{0} \right)}{\phi^{''}\left( r_{0} \right)}}\end{Bmatrix}} \right\rbrack}}}} \\{= {\int_{- \infty}^{\infty}\ {{\mathbb{d}u_{j}}{\exp\left\lbrack {{- {{\beta\phi}^{''}\left( r_{0} \right)}}\begin{Bmatrix}{\frac{\left( {u_{i}^{2} + u_{k}^{2}} \right)}{2} +} \\\begin{matrix}{u_{j}^{2} - {u_{j} \cdot}} \\{\left( {u_{j} + u_{k}} \right) +}\end{matrix} \\{2\frac{\phi\left( r_{0} \right)}{\phi^{''}\left( r_{0} \right)}}\end{Bmatrix}} \right\rbrack}}}} \\{= {\left( \frac{2\;\pi}{\beta\;{\phi^{''}\left( r_{0} \right)}} \right)^{\frac{3}{2}}{\exp\left\lbrack {{- 2}\;{{\beta\phi}^{''}\left( r_{0} \right)}\begin{Bmatrix}{\frac{\left( {u_{i} - u_{k}} \right)^{2}}{4.2} +} \\\frac{\phi\left( r_{0} \right)}{\phi^{''}\left( r_{0} \right)}\end{Bmatrix}} \right\rbrack}}}\end{matrix} & (11)\end{matrix}$

With the coarse graining, a lattice constant becomes α (=2) times, inthe same form as the original potential function on the basis of avariable u′ subjected to scale conversion by the equation u′=u/α. Thatis, the following similarity rule is obtained.

$\begin{matrix}\begin{matrix}{{\exp\left( {- {{\beta\phi}^{\prime}\left( {q_{i} - q_{k}} \right)}} \right)} = {\left( \frac{2\;\pi}{\beta\;{\phi^{''}\left( r_{0} \right)}} \right)^{\frac{3}{2}}{\exp\left\lbrack {{- {{\alpha\beta\phi}^{''}\left( r_{0} \right)}}\begin{Bmatrix}{\frac{\left( {u_{i}^{\prime} - u_{k}^{\prime}} \right)^{2}}{2} +} \\\frac{\phi\left( r_{0} \right)}{\phi^{''}\left( r_{0} \right)}\end{Bmatrix}} \right\rbrack}}} \\{= {\left( \frac{2\;\pi}{\beta\;{\phi^{''}\left( r_{0} \right)}} \right)^{\frac{3}{2}}{\exp\left( {{- {\alpha\beta ɛ}}\;{f\left( \frac{{q_{i}/\alpha} - {q_{k}/\alpha}}{\sigma} \right)}} \right)}}}\end{matrix} & (12)\end{matrix}$

However, in moving to the second row of the equation (12), it wasconfirmed that the same renormalization conversion is also establishedother than near the saddle point, α is called a first renormalizationfactor.

Coarse graining in a d-dimensional simple cubic lattice can be realizedby a potential moving method. Potential Moving is explained in, forexample, Leo P. Kadanoff, “STATISTICAL PHYSICS, Statics, Dynamics andRenormalization”, World Scientific (1999), Chap. 14.

A way of thinking potential moving will be described with reference toFIG. 1B. In potential moving, a multidimensional lattice is resolved toa one-dimensional chain. FIG. 1B shows an example of a two-dimensionallattice.

As shown in the upper part of FIG. 1B, particles are arranged on thelattice points of a two-dimensional cubic lattice. The interaction(nearest neighbor interaction) between nearest neighbor particles isindicated by a solid line. The direction parallel to the side of thelattice is referred to as the X direction, and the directionperpendicular to the X direction is referred to as the Y direction.

As shown in the middle part of FIG. 1B, it is considered that, for theparticles arranged in the X direction, every other particle isintegrated, i.e., a particle on a lattice point (lattice pointforcermation) indicated by an unshaded circle is integrated. A particleto be integrated (erased) is called an integrated particle. The nearestneighbor interaction between integrated particles is indicated by abroken line. If there is no nearest neighbor interaction betweenintegrated particles, the integrated particles are regarded as aone-dimensional chain extending in the X direction.

In potential moving, as shown in the lower part of FIG. 1B, the nearestneighbor interaction between the integrated particles is distributed tothe particles arranged on both side of an integrated particle in the Xdirection. The inter-particle interaction (double coupling) to which thedistributed interaction is added is indicated by a double line. Theinteraction indicated by the double line has intensity two times greaterthan the original nearest neighbor interaction. In this way, for theintegrated particles, a two-dimensional lattice can be converted to aone-dimensional chain. Thus, for the integrated particles, a method ofcoarse graining of the above-described one-dimensional chain can beexecuted. For example, in the case of a three-dimensional lattice, thesame procedure as described above is repeated in the Y and Z directions.In this way, coarse graining in a multidimensional lattice is executed.

If the potential moving method is used, for coarse graining of aninteraction portion of the partition function Z(β) in the d-dimensionalsimple cubic lattice, the following result is derived.

$\begin{matrix}{{\int{D_{q}^{N}{\exp\left( {{- \beta}{\sum\limits_{j = 1}^{N}{\sum\limits_{i = {j + 1}}^{N}{\phi\left( {q_{i} - q_{j}} \right)}}}} \right)}}} \propto {\int{D_{q}^{N^{\prime}}{\exp\left( {{- {\beta\alpha}^{d}}{\sum\limits_{j = 1}^{N}{\sum\limits_{i = {j + 1}}^{N}{ɛ\;{f\left( \frac{{q_{i}/\alpha} - {q_{j}/\alpha}}{\sigma} \right)}}}}} \right)}}}} & (13)\end{matrix}$

Here, the number N′ of renormalized particles and the firstrenormalization factor α are expressed by the following equations.N′=N/α ^(d)α=2^(n)

n is a positive integer and represents the number of times of coarsegraining.

Next, coarse graining of a kinetic energy portion of the partitionfunction Z(β) will be described. First, consideration is made in aone-dimensional chain. Particles j−1 and j which form a pair of nearestneighbor particles are coarsely grained to one particle. The way ofthinking of the block spin method of Kadanoff is applied to a momentum,and the momentum g_(j) of the coarsely grained particle having theparticles j and j−1 is defined as follows.g _(j)=(p _(j) +p _(j−1))/2

A weighting function T{P,p} is defined as follows.

${T\left\{ {P,p} \right\}} = {\prod\limits_{j = 1}^{\frac{N}{2}}\;{\delta\left( {P_{j} - g_{j}} \right)}}$

Here, δ is a δ-function of Dirac. Undoubtedly, the following conditionis satisfied.

${\int_{- \infty}^{\infty}{\prod\limits_{j = 1}^{\frac{N}{2}}\mspace{7mu}{{\mathbb{d}P_{j}}T\left\{ {P,p} \right\}}}} = 1$

Under this condition, the weighting function T{P,p} can be inserted intothe kinetic energy portion Z_(p) of the partition function, thus thefollowing equation is obtained.

$\begin{matrix}\begin{matrix}{Z_{p} = {\int{D_{p}^{\frac{N}{2}}{\int{D_{p}^{N}T\left\{ {P,p} \right\}{\exp\left( {{- \beta}{\sum\limits_{i}^{N}\frac{p_{i}^{2}}{2\; m}}} \right)}}}}}} \\{= {\int{D_{p}^{\frac{N}{2}}{\int{D_{p}^{N}{\exp\left( {{- \beta}{\sum\limits_{i}^{N}\frac{p_{i}^{2}}{2\; m}}} \right)}}}}}} \\{\prod\limits_{j = 1}^{\frac{N}{2}}\;{\delta\left( {P_{j} - \frac{p_{2\; j} + p_{{2\; j} - 1}}{2}} \right)}}\end{matrix} & (14)\end{matrix}$

With the integral representation of the δ-function, the kinetic energyportion Z_(p) of the partition function is expressed as follows.

$\begin{matrix}\begin{matrix}{Z_{p} = {\left( \frac{1}{2\;\pi} \right)^{\frac{N}{2}}{\int_{\infty}^{- \infty}{D_{p}^{\frac{N}{2}}{\int_{\infty}^{- \infty}{D_{t}^{\frac{N}{2}}{\int D_{p}^{N}}}}}}}} \\{\exp\left\lbrack {{{- \beta}{\sum\limits_{i}^{N}\frac{p_{i}^{2}}{2\; m}}} + {{\mathbb{i}}{\sum\limits_{j}^{\frac{N}{2}}{t_{j} \cdot \left( {P_{j} - \frac{p_{2\; j} + p_{{2\; j} - 1}}{2}} \right)}}}} \right\rbrack} \\{= {\left( \frac{1}{2\;\pi} \right)^{\frac{N}{2}}{\int_{\infty}^{- \infty}{D_{p}^{\frac{N}{2}}{\int{D_{t}^{\frac{N}{2}}{\exp\left( {{- {\mathbb{i}}}{\sum\limits_{i}^{\frac{N}{2}}{P_{j} \cdot t_{j}}}} \right)} \times}}}}}} \\{\int_{\infty}^{- \infty}{D_{p}^{N}{\exp\left\lbrack {- {\sum\limits_{j}^{\frac{N}{2}}\left( {{\frac{\beta}{2\; m}p_{2\; j}^{2}} + {\frac{\mathbb{i}}{2}{t_{j} \cdot p_{2\; j}}}} \right)}} \right\rbrack}}} \\{\exp\left\lbrack {- {\sum\limits_{j}^{\frac{N}{2}}\left( {{\frac{\beta}{2\; m}p_{{2\; j} - 1}^{2}} + {\frac{\mathbb{i}}{2}{t_{j} \cdot p_{{2\; j} - 1}}}} \right)}} \right\rbrack}\end{matrix} & (15)\end{matrix}$

This integration can be easily executed if the following formula isused.

${\int{D_{t}{\exp\left\lbrack {{\frac{1}{2}\left\langle {t,{Kt}} \right\rangle} + \left\langle {p,t} \right\rangle} \right\rbrack}}} = {\frac{1}{\sqrt{\det\left( \frac{K}{2\;\pi} \right)}}\exp\frac{1}{2}\left\langle {p,{K^{- 1}p}} \right\rangle}$

Here, <A,B> represents the inner product of vectors A and B. If thisformula is used, the kinetic energy portion Z_(p) of the partitionfunction can be expressed as follows.

$\begin{matrix}\begin{matrix}{Z_{p} = {\left( \frac{1}{2\;\pi} \right)^{\frac{3N}{2}}\left( \frac{2\;\pi\; m}{\beta} \right)^{\frac{3\; N}{2}}{\int{D_{p}^{\frac{N}{2}}{\int_{\infty}^{- \infty}D_{t}^{\frac{N}{2}}}}}}} \\{\exp\left( {{{- \frac{1}{2}}{\sum\limits_{j}^{\frac{N}{2}}{{t_{j} \cdot \left( \frac{m}{\beta} \right)}t_{j}}}} + {\sum\limits_{j}^{\frac{N}{2}}{{\mathbb{i}}\;{t_{j} \cdot P_{j}}}}} \right)} \\{= {\left( \frac{2\;\pi\; m}{\beta} \right)^{\frac{3\; N}{4}}{\int{D_{p}^{\frac{N}{2}}{\exp\left( {{- 2}\;\beta{\sum\limits_{j = 1}^{\frac{N}{2}}\frac{P_{j}^{2}}{2\; m}}} \right)}}}}}\end{matrix} & (16)\end{matrix}$

With this coarse graining, the kinetic energy becomes α (=2) times. Ifthe entire kinetic energy of two particles is divided into the kineticenergy of the center of gravity and the kinetic energy of the relativemomentum, and the relative momentum is integrated, the same result isobtained.

Next, consideration is made as to coarse graining in a d-dimensionalsimple cubic lattice. First, in the case of coarse graining in thedirection of the crystal axis a, the kinetic energy E_(kin) of aparticle becomes α times. Successively, in the case of coarse grainingin the b-axis direction, the kinetic energy E_(kin)α of a particlebecomes α times as large as E_(kin)α, that is, E_(kin)α². Thus, in thed-dimensional lattice, the kinetic energy becomes E_(kin)α^(d), and thefollowing relation expression is obtained.

$\begin{matrix}{{\int{D_{p}^{N}{\exp\left( {{- \beta}{\sum\limits_{j = 1}^{N}\frac{p_{j}^{2}}{2\; m}}} \right)}}} \propto {\int{D_{p}^{N^{\prime}}{\exp\left( {{- \beta}\;\alpha^{d}{\sum\limits_{j = 1}^{N^{\prime}}\frac{p_{j}^{2}}{2\; m}}} \right)}}}} & (17)\end{matrix}$

In the case of coarse graining n times, it is obvious that α=2^(n).

From the above discussion, the partition function Z(β) of the entireHamiltonian subjected to coarse graining can be expressed, excluding acoefficient which does not affect the result.

$\begin{matrix}{{Z(\beta)} = {\int{{\mathbb{d}\Gamma_{N^{\prime}}}{\exp\left\lbrack {{- {\beta\alpha}^{d}}{\sum\limits_{j = 1}^{N^{\prime}}\left\{ {\frac{p_{j}^{2}}{2m} + {ɛ{\sum\limits_{i = {j + 1}}^{N^{\prime}}{f\left( \frac{{q_{i}/\alpha} - {q_{j}/\alpha}}{\sigma} \right)}}}} \right\}}} \right\rbrack}}}} & (18)\end{matrix}$

Next, in order to set the Hamiltonian of the same form for the systemsubjected to coarse graining, “renormalized” variables are introduced.As described below, a renormalized position vector q′, a renormalizedmomentum vector p′, a renormalized mass m′, a renormalized interactioncoefficient ε′, and a renormalized coefficient β′ are defined. Scaleconversion is done such that the form of the Hamiltonian is invariable.q′=q/αp′=pα ^(δ/2)m′=mα ^(δ−γ)ε′=εα^(γ)β′=βα^(d−γ)

Here, γ is equal to or greater than 0 and equal to or smaller than d. γis called a second renormalization factor. δ has a value equal to orgreater than 0. δ is called a third renormalization factor.

The renormalized variables are defined in the above-described manner,such that the renormalized Hamiltonian H′ of the same form as theHamiltonian H is obtained. The renormalized Hamiltonian H′ is expressedas follows.

$\begin{matrix}{H^{\prime} = {\sum\limits_{j = 1}^{N^{\prime}}\left\{ {\frac{p_{j}^{\prime 2}}{2m^{\prime}} + {ɛ^{\prime}{\sum\limits_{i = {j + 1}}^{N^{\prime}}{f\left( \frac{q_{i}^{\prime} - q_{j}^{\prime}}{\sigma} \right)}}}} \right\}}} & (19)\end{matrix}$

From the equation (18), the partition function is as follows.Z(β)=∫dΓ _(N′)′exp(−β′H′(p′,q′;m′,e′))  (20)

Here, a volume element of a phase space is as follows, excluding acoefficient which does not affect the result.

$\begin{matrix}{{\mathbb{d}\Gamma_{N^{\prime}}^{\prime}} = {\frac{1}{W_{N^{\prime}}}{\prod\limits_{j = 1}^{N^{\prime}}\;{{\mathbb{d}p_{j}^{\prime}}{\mathbb{d}q_{j}^{\prime}}}}}} & (21)\end{matrix}$

In the system (renormalized particle system S′) which is described bythe renormalized Hamiltonian H′, the equations of motion are given fromthe following canonical equations.

$\begin{matrix}{\frac{\mathbb{d}p_{i}^{\prime}}{\mathbb{d}t^{\prime}} = {\frac{\partial H^{\prime}}{\partial q_{i}^{\prime}}\mspace{14mu}\left( {{i = 1},\ldots\mspace{14mu},N^{\prime}} \right)}} & (22) \\{\frac{\mathbb{d}q_{i}^{\prime}}{\mathbb{d}t^{\prime}} = {\frac{\partial H^{\prime}}{\partial p_{i}^{\prime}}\mspace{14mu}\left( {{i = 1},\ldots\mspace{14mu},N^{\prime}} \right)}} & (23)\end{matrix}$

The equations of motion become the following equations.

$\begin{matrix}{\frac{\mathbb{d}p_{i}^{\prime}}{\mathbb{d}t^{\prime}} = {{- ɛ^{\prime}}{\sum\limits_{j \neq i}^{N^{\prime}}{\frac{\partial}{\partial q_{i}^{\prime}}{f\left( \frac{q_{i}^{\prime} - q_{j}^{\prime}}{\sigma} \right)}\mspace{14mu}\left( {{i = 1},\ldots\mspace{14mu},N} \right)}}}} & (24) \\{\frac{\mathbb{d}q_{i}^{\prime}}{\mathbb{d}t^{\prime}} = {\frac{p_{i}^{\prime}}{m^{\prime}} = {v_{i}^{\prime}\mspace{14mu}\left( {{i = 1},\ldots\mspace{14mu},N^{\prime}} \right)}}} & (25)\end{matrix}$

Here, the following equation is established such that a renormalizedtime t′ matches with the renormalized variables q′, p′, m′, and ε′.t′=q′/v′t/α ^((1+γ−δ/2))

Molecular dynamics calculation can be carried out on the renormalizedparticle system S on the basis of the equations (24) and (25) of motion.

The above consideration is concluded. The particle system S and therenormalized particle system S′ are associated with each other by ascale conversion rule expressed by the following conversion equations.

The number N of particles of the particle system S and the number N′ ofrenormalized particles of the renormalized particle system S′ areassociated with each other by the following conversion equation.N′=N/α ^(d)  (conversion equation N)

The position vector q of a certain particle of the particle system S andthe position vector q′ of a certain particle of the renormalizedparticle system S′ are associated with each other by the followingconversion equation.q′=q/α  (conversion equation q)

The momentum vector p of a certain particle of the particle system S andthe momentum vector p′ of a certain particle of the renormalizedparticle system S′ are associated with each other by the followingconversion equation.p′=pα ^(δ/2)  (conversion equation p)

The mass m of each particle of the particle system S and the mass m′ ofeach particle of the renormalized particle system S are associated witheach other by the following conversion equation.m′=mα ^(δ−γ)  (conversion equation m)

The interaction coefficient ε in the particle system S and theinteraction coefficient ε′ in the renormalized particle system S′ areassociated with each other by the following conversion equation.ε′=εα^(γ)  (conversion equation ε)

A time scale t (a certain time interval t in molecular dynamicscalculation which is supposed to be carried out on the particle systemS) when the equation of motion is resolved for the particle system S anda time scale t′ (a certain time interval t′ in molecular dynamicscalculation which is carried out on the renormalized particle system S′)when the equation of motion is resolved for the renormalized particlesystem S′ are associated with each other by the following conversionequation.t′=t/α ^((1+γ−δ/2))  (conversion equation t)

The velocity vector v of a certain particle of the particle system S andthe velocity vector v′ of a certain particle of the renormalizedparticle system S′ are associated with each other by the followingconversion equation because the velocity vector v is expressed by v=p/m.v′=vα ^((γ−δ/2))  (conversion equation v)

Here, the dimension number of a space where the particle system S isarranged is d. The first renormalization factor α is 2^(n) (where n is apositive integer), and the second renormalization factor γ is equal toor greater than 0 and equal to or smaller than d. The thirdrenormalization factor δ is equal to or greater than 0.

When the value of the second renormalization factor γ is large, forexample, equal to d, a very large renormalized interaction coefficientε′ and a small renormalized mass m′ appear in numerical calculation. Innumerical calculation, it is preferable that the second renormalizationfactor γ is set to 0. At this time, if the third renormalization factorδ is also set to 0, only the time t is converted, and the MD of therelated art can be used.

For example, when δ=2, with regard to a spring-mass system of aone-dimensional chain having a mass m and a spring constant s, thedispersion relationship is obtained as follows from the renormalizedHamiltonian.

$\begin{matrix}{\omega = {{\pm \left( \frac{4\; s}{m} \right)^{\frac{1}{2}}}\frac{1}{\alpha}{\sin\left( \frac{k\;\alpha\; a}{2} \right)}}} & (26)\end{matrix}$

Here, an angular velocity is ω, and a wave number is k. An inter-atomdistance in an equilibrium state is a. The second renormalization factorγ was set to 0. The equation (26) becomes as follows in the limit of along wavelength, and has no relation to the first renormalization factorα.

$\begin{matrix}{\omega = {{\pm \left( \frac{{sa}^{2}}{m} \right)^{\frac{1}{2}}}k}} & (27)\end{matrix}$

Next, an example of a computer simulation by renormalization groupmolecular dynamics will be described.

First, a simulation in which heat flows in aluminum to raise thetemperature, and a melting point and latent heat are calculated will bedescribed. In a first simulation, for comparison, molecular dynamicscalculation (with no renormalization applied) of the related art wascarried out on a system having 16000 aluminum atoms. The aluminum atomsare arranged in a three-dimensional space, and the dimension number d is3.

In second and third simulations, molecular dynamics calculation wascarried out on a system renormalized from the system corresponding tothe first simulation. The first renormalization factor α is determinedto be 2 in the second simulation and 2² (=4) in the third simulation.From the (conversion equation N), the number of particles of the systemin the second simulation becomes 2000, and the number of particles ofthe system in the third simulation becomes 250. The secondrenormalization factor γ was determined to be 0. The thirdrenormalization factor δ was determined to be 2. The related art method(first simulation) corresponds to a case where the first renormalizationfactor α is set to 1.

In the first simulation, the Morse potential was used as an inter-atompotential (that is, an inter-particle potential in the original system).This potential is expressed as follows.φ(r)=ε{e ^(−2A(r−r) ⁰ ⁾−2e ^(−A(r−r) ⁰ ^()})  (28)

Here, the interaction coefficient ε is 1.92×10⁻²⁰ [J], the constant A is2.35×10¹⁰ [m⁻¹], and the constant r₀ is 2.86×10⁻¹⁰ [m]. Theinter-particle distance is r.

In the second and third simulations, the same form as in the equation(28) is used as the inter-particle potential (that is, theinter-particle potential in the renormalized system). The renormalizedinteraction coefficient in each of the second and third simulations isfound by the (conversion equation ε). Since the second renormalizationfactor γ is 0, the interaction coefficient in each of the second andthird simulations is the same as the interaction coefficient εin thefirst simulation. The constant r₀ and the constant A are the same asthose in the first simulation.

The mass of each aluminum atom in the first simulation (that is, themass of each particle in the original system) is 4.48×10⁻²⁶ [kg]. Therenormalized mass of each particle in each of the second and thirdsimulations is found by the (conversion equation m). The mass in thesecond simulation becomes 17.9×10⁻²⁶ [kg] which is obtained bymultiplying the mass in the first simulation by 2² (=4), and the mass inthe third simulation becomes 71.7×10⁻²⁶ [kg] which is obtained bymultiplying the mass in the first simulation by 4² (=16).

Next, the initial arrangement of particles will be described. In theinitial state, the particles are arranged on the lattice points of aface-centered cubic lattice. In the first simulation, 40 layers, in eachof which 20 particles are arranged in a width direction and 20 particlesare arranged in a depth direction, are arranged to overlap each other. Aportion of 1/α in the width direction, 1/α in the depth direction, and1/α in the height direction of the original system (a portion of 1/α³ ofthe original system) corresponds to the renormalized system. In thesecond simulation, 20 layers in each of which 10 particles are arrangedin the width direction and 10 particles are arranged in the depthdirection are arranged to overlap each other. In the third simulation,10 layers in each of which 5 particles are arranged in the widthdirection and 5 particles are arranged in the depth direction arearranged to overlap each other.

Next, a time width of numerical integration will be described. It isassumed that a time width in the original system is Δt, and the velocityof the particle is v. The time width Δt is determined, for example, soas to satisfy the following equation (“<<” means that the left side ofthe equation is sufficiently smaller than the right side).

$\frac{v\;\Delta\; t}{r_{0}}{\operatorname{<<}1}$

Here, r₀ is a constant in the inter-particle potential, and gives anindication of the inter-particle distance. It is assumed that a timewidth in the renormalized system is Δt′, and the velocity of theparticle is v′. From the (conversion equation v), since v′=v/α, the timewidth Δt′ in the renormalized system is set to, for example, αΔt.

The time width of numerical integration was set to 5.0 [fs] in the firstsimulation, 10 [fs] in the second simulation, and 20 [fs] in the thirdsimulation.

Next, an initial velocity which is given to a particle will bedescribed. First, the relationship between the temperature of theoriginal system and the temperature of the renormalized system will bedescribed. The temperature T of the original system is defined by thefollowing equation.

$\begin{matrix}{{\frac{3}{2}k_{B}T} = {\frac{1}{N - 1}{\sum\limits_{i}^{N}{\frac{1}{2}{m\left( {v_{i} - v_{g}} \right)}^{2}}}}} & (29)\end{matrix}$

Here, k_(B) is the Boltzmann constant, v_(i) is the velocity vector of aparticle i, and v_(g) is the velocity vector of the center of gravity ofthe particle system. The temperature is, that is, the “fluctuation” ofthe velocity. The temperature T′ of the renormalized system is definedby the following equation.

$\begin{matrix}{{\frac{3}{2}k_{B}T^{\prime}} = {\frac{1}{N^{\prime} - 1}{\sum\limits_{i}^{N^{\prime}}{\frac{1}{2}{m^{\prime}\left( {v_{i}^{\prime} - v_{g}^{\prime}} \right)}^{2}}}}} & (30)\end{matrix}$

In the original system, the kinetic energy of a particle having a mass mand a velocity v is ½ mv². In the renormalized system, the kineticenergy of a particle having a mass m′ and a velocity v′ is ½ m′v′². Fromthe (conversion equation m) and the (conversion equation v), thefollowing relationship is obtained.

$\begin{matrix}{{\frac{1}{2}{mv}^{2}} = {{\frac{1}{2}\frac{m}{\alpha^{2}}\left( {v\;\alpha} \right)^{2}} = {\frac{1}{2}m^{\prime}v^{\prime\; 2}}}} & (31)\end{matrix}$

Thus, the temperature T of the original system becomes equal to thetemperature T′ of the renormalized system. In the first to thirdsimulations, the initial velocity is given such that the temperature ofthe particle system is at 300 [K].

Next, a method which flows heat into the particle system will bedescribed. In the first to third simulations, the kinetic energy wasgiven to particles arranged at the bottom of the particle system, andthe inflow of heat into the system was simulated. A surface into whichheats flows is called a heated surface.

First, a heat inflow method in the first simulation (that is, a heatinflow method in the original system) will be described. It is assumedthat kinetic energy ΔE is given per particle on the heated surface witha time width Δt. A variation in the momentum vector of the particle isobtained by the following equation (32). Here, the momentum vector at acertain time t is p_(t), and the momentum vector at a time t+Δt isP_(t+Δt).

$\begin{matrix}{{\frac{p_{t}^{2}}{2\; m} - \frac{p_{t + {\Delta\; t}}^{2}}{2\; m}} = {\Delta\; E}} & (32)\end{matrix}$

It is assumed that the momentum vector becomes λ times because of thekinetic energy ΔE. That is, the following equation is assumed.p _(t+Δt) =λp _(t)  (33)

If the equation (33) is substituted in the equation (32), λ is expressedas follows.

$\begin{matrix}{\lambda = \sqrt{1 + \frac{2\; m\;\Delta\; E}{p_{t}^{2}}}} & (34)\end{matrix}$

The momentum vector is updated by the equation (33) for every time widthΔt. Thus, the inflow of heat can be simulated such that the kineticenergy ΔE flows per time width Δt.

If the number of particles on the heated surface is N_(h), and the areaof the heated surface is S_(h), power f_(h) which flows per unit area isexpressed as follows.

$\begin{matrix}{f_{h} = \frac{N_{h}\Delta\; E}{S_{h}\Delta\; t}} & (35)\end{matrix}$

In the first simulation, the number N_(h) of particles on the heatedsurface is 400, the area S_(h) of the heated surface is 2.83×10⁻¹⁷ [m²],the kinetic energy ΔE is 4.14×10⁻²⁴ [J], the time width Δt is 5.0 [fs],and the power f_(h) is 1.17×10¹⁰ [W/m²].

In the second and third simulations (in the renormalized system), powerf_(h)′ which flows per unit area is given by the following equation.f _(h) ′=f _(h)/α  (36)

The coefficient λ′ for updating the momentum vector is given by thefollowing equation.

$\begin{matrix}{\lambda^{\prime} = \sqrt{1 + \frac{2\; m^{\prime}\Delta\; E}{p_{t^{\prime}}^{\prime\; 2}}}} & (37)\end{matrix}$

From the equation (31), the kinetic energy ΔE to be given becomes equalbetween the original system and the renormalized system.

Next, simulation results will be described with reference to FIG. 2. Ofthe graphs shown in FIG. 2, a lower graph shows the result of the firstsimulation (molecular dynamics calculation of the related art), a middlegraph shows the result of the second simulation (renormalization groupmolecular dynamics calculation when the renormalization factor α is 2),and an upper graph shows the result of the third simulation(renormalization group molecular dynamics calculation when therenormalization factor α is 4). In the three graphs, the vertical axisrepresents temperature (unit [K]), and the horizontal axis representsthe amount of inflow heat, that is, inflow energy (unit 10⁻¹⁶ [J]).

First, the result of the first simulation will be described withreference to the lower graph. When no energy flows into the system, thetemperature of the system is 300 [K]. As the inflow energy increases,the temperature of the system rises. If the inflow energy reaches about5×10⁻¹⁶ [J], temperature rise temporarily stops at about 850 [K], andthe temperature is substantially constant until the inflow energyreaches about 9×10⁻¹⁶ [J]. This temperature represents a melting point.If the inflow energy exceeds about 9×10⁻¹⁶ [J], the temperature risesagain.

Energy which flows into the system during a period in which thetemperature is constant corresponds to latent heat. Latent heat which isobtained from this simulation becomes 14.7 [kJ/mol]. The melting pointof aluminum experimentally obtained is 933 [K], and latent heat is 11[kJ/mol]. It is considered that the difference between the valueobtained by the simulation and an experimental value is due to theprecision of the used inter-atom potential or the like.

Next, the results of the second and third embodiments will be describedwith reference to the middle and upper graphs. The results of the secondand third simulations have the same tendency as the result of the firstsimulation. That is, there is a tendency that the temperature issubstantially constant around at 850 [K] during a period in which theinflow energy is about 5×10⁻¹⁶ [J] to about 9×10⁻¹⁶ [J]. Thus, from thesecond and third simulations, it can be obtained that the melting pointis about 850 [K], and latent heat is 14.7 [kJ/mol].

The energy which flows into the particle system is obtained byintegrating the inflow power by the period in which the simulation isdone. If the energy which flows into the renormalized system is E′, theenergy E which flows into the original system is obtained by thefollowing equation.E=α ³ E′  (38)

The inflow energy per particle is equal between the original system andthe renormalized system.

The temperature of the renormalized system is obtained from the equation(30). From the above discussion, the temperature of the original systemis equal to the temperature of the renormalized system.

Although in the above-described simulation, the inflow energy andtemperature in the original system are obtained, if necessary, otherphysical quantities in the original system, for example, the position,momentum, velocity, and the like of the particle can be obtained on thebasis of the above-described scale conversion rule.

As described above, in renormalization group molecular dynamics, thephysical quantities (in the above-described example, the melting pointand latent heat) which are expected to be obtained by molecular dynamicscalculation on a system before renormalization (the system correspondingto the first simulation in the above-described example) can becalculated by carrying out molecular dynamics calculation on a systemhaving a smaller number of particles than the system beforerenormalization. Various kinds of calculation necessary for a simulationby renormalization group molecular dynamics can be carried out on acomputer on the basis of a program.

In the simulation of the example, aluminum has a face-centered cubiclattice. Although the scale conversion rule is derived from theconsideration on a particle system arranged on a simple cubit lattice,the invention is not limited to a simple cubic structure and iseffectively applied to other crystal structures.

Next, a simulation in which heat flows into silicon to raise thetemperature, and a melting point and latent heat are derived will bedescribed. Silicon has a diamond structure. Similarly to the simulationof aluminum, three simulations were carried out. A first simulation is amolecular dynamics method of the related art for comparison. In secondand third simulations, molecular dynamics calculation was carried out ona system renormalized from the system corresponding to the firstsimulation. The renormalization factor α was set to 2 in the secondsimulation and 2² (=4) in the third simulation. The number of particlesis 8192 in the first simulation, 1024 in the second simulation, and 128in the third simulation.

As the inter-atom potential in the first simulation (that is, theinter-particle potential in the original system), the S-W(Stillinger-Weber) potential was used. This potential is expressed.

$\begin{matrix}{\phi = {{\sum\limits_{i < j}{ɛ\;{f\left( {r_{ij}/\sigma} \right)}}} + {\sum\limits_{i < j < k}{ɛ\;{g\left( {{r_{i}/\sigma},{r_{j}/\sigma},{r_{k}/\sigma}} \right)}}}}} & (39)\end{matrix}$

Here, the definition is made as follows.f(x)=A(B/x ^(P)−1/x ^(q))exp{1/(x−a)} (x<a)f(x)=0 (x>a)g(x _(i) ,x _(j) ,x _(k))=h(x _(ij) ,x _(ik),θ_(ijk))+h(x _(ji) ,x_(jk),θ_(jik))+h(x _(ki) ,x _(kj) ,θ _(kij))h(x _(ij) ,x _(ik),θ_(ijk))=λ exp{γ/(x _(ij) −a)+γ/(x _(ik)−a)}(cosθ_(ijk)+⅓)² (x _(ij) <a,x _(ik) <a)h(x _(ij) ,x _(ik),θ_(ijk))=0 (x _(ij) >a,x _(ik) >a)

The interaction coefficient ε is 2.167 [eV], σ is 2.095 [Å], theconstants A, B, p, q, a, λ, and γ are respectively 7.05, 0.602, 4.0,0.0, 1.8, 21.0, and 1.2.

As shown in FIG. 3A, the angle of coupling between particles i and k andcoupling between particles i and j is θ_(ijk), and the inter-particledistance from the particle i to the particle j is r_(ij). Therenormalized interaction coefficient in the inter-particle potential ofeach of the second and third simulations is obtained by the same methodas in the simulation of aluminum.

The mass of a silicon atom in the first simulation (that is, the mass ofa particle in the original system) was set to 4.67×10⁻²⁶ [kg]. Therenormalized mass in each of the second and third simulations isobtained by the same method as in the simulation of aluminum.

Next, simulation results will be described with reference to FIG. 3B. Ofthe graphs shown in FIG. 3B, a lower graph shows the result of the firstsimulation (molecular dynamics calculation of the related art), a middlegraphs shows the result of the second simulation (renormalization groupmolecular dynamics calculation when the renormalization factor α is 2),and an upper graph shows the result of the third simulation(renormalization group molecular dynamics calculation when therenormalization factor α is 4).

First, the result of the first simulation will be described withreference to the lower graph. During a period in which the inflow energyincreases from about 1.3×10⁻¹⁶ [J] to about 4.6×10⁻¹⁶ [J], thetemperature is substantially constant at about 1750 [K]. Latent heat isobtained to be 25.6 [kJ/mol]. The melting point of siliconexperimentally obtained is 1687 [K], and latent heat is 50.6 [kJ/mol].It is considered that the difference between the value obtained in thesimulation and the experimental value is due to the precision of theinter-atom potential or the like.

Next, the results of the second and third simulations will be describedwith reference to the middle and upper graphs. The results of the secondand third simulations substantially have the same tendency as the resultof the first simulation. Thus, the melting point and latent heat closeto the values obtained in the first simulation can be obtained on thebasis of the second and third simulations.

Next, a simulation in which a dispersion relationship is obtained byrenormalization group molecular dynamics will be described withreference to FIG. 4. Vibration was applied in parallel to the [111]direction of aluminum to calculate the wave number in the [111]direction. The atoms were connected to each other through springs. Aspring constant s was set to 40.3 [N/m]. The spring constant s isobtained by approximating the near-minimum value of the inter-atompotential of aluminum using a quadratic function. The mass m of analuminum atom was set to 4.48×10⁻²⁶ [kg]. The inter-atom distance a wasset to the 2.86×10⁻¹⁰ [m]. The exact solution (the solution analyticallyobtained) of the dispersion relationship is given by the followingequation when an oscillatable angular frequency is ω and the wave numberin the [111] direction is k.

$\omega = {\left( \frac{4s}{m} \right)^{1/2}{\sin\left( \frac{ka}{2} \right)}}$

The exact solution can be obtained by α=1 in the equation (26). A solidline L in the graph indicates an exact solution. A point P1 indicatesthe result of molecular dynamics calculation of the related art(corresponding to when the renormalization factor α=1). Points P2 to P5respectively indicate the results when the renormalization factor α is2⁷, 2¹⁵, 2²¹, and 2²⁸. The point that the renormalization factor αincreases means that the size of the original system (the system beforerenormalization) increases. The renormalization factor α increases,obtaining a small wave number. The result obtained by renormalizationgroup molecular dynamics satisfactorily coincides with an exactsolution. In this way, if renormalization group molecular dynamics isused, calculation in a long wavelength region which is not easilycarried out in molecular dynamics of the related art (calculation on amacro system which is regarded as a continuous body) can be carried outwith good precision.

An XYZ orthogonal coordinate system is considered in a space where aparticle system is arranged. In the above-described example, therenormalized position vector q′ was obtained from the position vector qof the original system on the basis of the renormalization factor α.q′=q/α  (conversion equation q)

This indicates that spatial scale conversion in all the X, Y, and Zdirections is carried out by the common renormalization factor α. Therenormalization factor α can be set separately in each of the X, Y, andZ directions. That is, scale conversion having anisotropy may be carriedout.

When the second renormalization factor γ is determined to be 0 and thethird renormalization factor α is determined to be 2, a scale conversionrule having anisotropy is expressed by the following conversionequation. The renormalization factor in the X direction is α_(X), therenormalization factor in the Y direction is α_(Y), and therenormalization factor in the Z direction is α_(Z). Each of therenormalization factors α_(X) to α_(Z) is, for example, 2^(n) (where nis a positive integer).

The number N of particles of the particle system S and the number N′ ofrenormalized particles of the renormalized particle system S′ areassociated with each other by the following conversion equation.N′=N/α _(X)α_(Y)α_(Z)  (conversion equation Na)

Each component of the position vector q of a certain particle of theparticle system S and each component of the position vector q′ of acertain particle of the renormalized particle system S′ are associatedwith each other by the following conversion equation.q′ _(η) =q _(η)/α_(η)  (conversion equation qa)

Here, η represents one of the X, Y, and Z components of a vector(hereinafter, the same is applied to (conversion equation pa),(conversion equation va), and (conversion equation ma)). Each componentof the momentum vector p of a certain of the particle system S and eachcomponent of the momentum vector p′ of a certain particle of therenormalized particle system S′ are associated with each other by thefollowing conversion equation.p _(η) ′=p _(η)α_(η)  (conversion equation pa)

When γ=0 and δ is undetermined, δ is determined, and the followingequation is substituted, such that association is made.p _(η) ′=p _(η)α_(η) ^(δ/2)

Each component of the velocity vector v of a certain particle of theparticle system S and each component of the velocity vector v′ of acertain particle of the renormalized particle system S are associatedwith each other by the following conversion equation.v _(η) ′=v _(η)/α_(η)  (conversion equation va)

When γ=0 and δ is undetermined, δ is determined, and the followingequation is substituted, such that association is made.v _(η) ′=v _(η)/α_(η) ^(δ/2)

Since the renormalization factor in each direction is set, the mass ofthe particle is subjected to scale conversion in each direction. Themass m of each particle of the particle system S and the mass m′_(X) inthe X direction, the mass m′_(Y) in the Y direction, and the mass m′_(Z)in the Z direction of each particle of the renormalized particle systemS′ are associated with each other by the following conversion equation.m _(η) ′=mα _(η) ²  (conversion equation ma)

The corresponding mass is used by each direction in the equation (25) ofmotion in the renormalized system.

When γ=0 and δ is undetermined, δ is determined, and the followingequation is substituted, such that association is made.m _(η) ′=mα _(η) ^(δ)

The interaction coefficient ε in the particle system S and theinteraction coefficient ε′ in the renormalized particle system S′ areassociated with each other by the following conversion equation.ε′=ε  (conversion equation εa)

The time scale t in the particle system S and the scale t′ in therenormalized particle system S are associated with each other by thefollowing conversion equation.t′=t  (conversion equation ta)

For example, a simulation of a thin film is considered. A thin film hasa larger film thickness in a direction parallel to the film surface thanin a film thickness direction. Thus, the renormalization factor isgreater in the direction parallel to the film surface than in the filmthickness direction. For example, if the film thickness direction is theZ direction, it should suffice that the renormalization factors α_(X)and α_(Y) in the direction parallel to the film surface are set to begreater than the renormalization factor α_(Z) in the film thicknessdirection.

The second scaling parameter γ is included in the (conversion equationm), the (conversion equation ε), and the (conversion equation t) in theform of α^(γ). When γ=0 is set, α^(γ) simply becomes “1”. Thus, forexample, when it is determined to be δ=2, with regard to the mass of theparticle, the interaction coefficient, and the time, the following scaleconversion may be carried out without using the second scaling parameterγ (without introducing a second scaling parameter as γ).m′=mα ²ε′=εt′=t

In this case, the same scale conversion as when the second scalingparameter γ is set to 0 is carried out. When γ=0 and δ=2, with regard tothe velocity vector, the following scale conversion is carried out.v′=v/α

If the renormalization factors α, α_(X), α_(Y), and α_(Z) are greaterthan 1, the number of particles in the renormalized system is smallerthan the number of particles in the original system.

If necessary, the renormalization factor (for example, α) may differbetween a certain portion and another portion of a system which will besubjected to a simulation.

In the above-described example, as the example of renormalization groupmolecular dynamics, the melting point, latent heat, and the dispersionrelationship have been obtained. Renormalization group moleculardynamics is not limited thereto and may be used for the simulations ofvarious kinds of physical phenomenon. Renormalization group moleculardynamics uses molecular dynamics and can be thus applied to variousproblems which are handled in molecular dynamics. For example, problemsof flow, structure, heat, reaction, and the like, problems of mechanism,abrasion, and lubrication, and the like can be handled. Withrenormalization group molecular dynamics, it is possible to obtain thephysical quantities of a system having a larger number of particles thana system in which molecular dynamics calculation is actually carriedout. Thus, if renormalization group molecular dynamics which has beenfound by the inventors is used, it becomes possible to carry out aprecise simulation on a problem of scaling on which a simulation is noteasily carried out in the related art. As the method of moleculardynamics calculation which is executed in renormalization groupmolecular dynamics, various known methods may be used.

In NPL 1 described in the “DESCRIPTION OF THE RELATED ART” of thisspecification, there is an attempt to apply the renormalizationtechnique to molecular dynamics. However, in NPL 1, only the interactionportion of the Hamiltonian is renormalized, and the kinetic energyportion is not renormalized. In NPL 1, with regard to the momentum (orvelocity), simple average, not renormalization, is taken. For thisreason, static analysis is still done, a temperature may not be takeninto consideration explicitly, and an artificial method (the equations(23), (24), and (25) described in “3. Numerical Calculation Method” ofNPL 1) is necessary for forcibly causing energy dissipation. For thisreason, at the time of the application to calculation of finitetemperature or dynamic problems involved in vibration, noise, and phasetransition, there is a problem in that reliability is lacking. The scaleconversion rule described in NPL 1 corresponds to a case when γ=0 andδ=3 in the invention.

In the course of deriving the scale conversion rule of the inventors,the kinetic energy portion of the Hamiltonian is renormalized. Thus, anappropriate scale conversion rule of the momentum vector (or velocityvector) is derived. Therefore, it becomes possible to take intoconsideration a temperature explicitly in a simulation.

Although the invention has been described in connection with theexample, the invention is not limited to the example. For example, it isobvious to those skilled in the art that various changes, improvements,combinations, or the like may be made.

What are claimed are:
 1. A simulation method comprising the steps of:(a) determining, via a processor, a first renormalization factor αgreater than 1, a second renormalization factor γ equal to or greaterthan 0 and equal to or smaller than d, and a third renormalizationfactor δ equal to or greater than 0 using a dimension number d of aspace where a particle system S is arranged to obtain a number N′ ofrenormalized particles by a conversion equation N′=N/α^(d), to obtain amass m′ of each renormalized particle by a conversion equationm′=mα^(δ)/α^(γ), and to obtain a renormalized interaction coefficient ε′by a conversion equation ε′=εα^(γ), wherein the particle system S inwhich a number of particles is N, a mass of each particle is m, andinter-particle interaction potential energy can be expressed by aproduct εf of a dimensionless function f representing an inter-particledistance dependence and an interaction coefficient ε; and (b)performing, via the processor, a molecular dynamics calculation on arenormalized particle system S′ in which the number of renormalizedparticles is N′, the mass of each particle is the mass m′ of eachrenormalized particle, and inter-particle interaction potential energyis expressed by a product ε′f of the dimensionless function f and therenormalized interaction coefficient ε′.
 2. The simulation methodaccording to claim 1, further comprising the step of: (c) when theposition vector of a certain particle in the renormalized particlesystem S′ obtained by molecular dynamics calculation in the step (b) isq′, the momentum vector of a certain particle in the renormalizedparticle system S′ obtained by molecular dynamics calculation in thestep (b) is p′, the velocity vector of a certain particle in therenormalized particle system S′ obtained by molecular dynamicscalculation in the step (b) is v′, and a certain time interval of acertain particle in the renormalized particle system S′ during moleculardynamics calculation in the step (b) is t′, carrying out at least one ofcalculation for obtaining the position vector q of a certain particle inthe particle system S by a conversion equation q=q′α, calculation forobtaining the momentum vector p of a certain particle in the particlesystem S by a conversion equation p=p′/α^(δ/2), calculation forobtaining the velocity vector v of a certain particle in the particlesystem S by a conversion equation v=v′α^((−γ+δ/2)), and calculation forobtaining a certain time interval t during molecular dynamicscalculation assumed to be carried out on the particle system S by aconversion equation t=t′^((1+γ−δ/2)) using at least one of the firstrenormalization factor α, the second renormalization factor γ, and thethird renormalization factor δ determined in the step (a).
 3. Thesimulation method according to claim 1, wherein the secondrenormalization factor γ is
 0. 4. The simulation method according toclaim 1, wherein the first renormalization factor is 2 ^(n) where n is apositive integer.
 5. The simulation method according to claim 1, whereinthe dimension number d is
 3. 6. A simulation method comprising the stepsof: (a) determining, via a processor, a first renormalization factor agreater than 1 and a third renormalization factor δ equal to or greaterthan 0 using a dimension number d of a space where a particle system Sis arranged to obtain the number N′ of renormalized particles byconversion equation N′=N/α^(d) and to obtain the mass m′ of eachrenormalized particle by a conversion equation m′=mα^(δ), wherein theparticle system S in which a number of particles is N, a mass of eachparticle is m and inter-particle interaction potential energy can beexpressed by a product εf of a dimensionless function f representinginter-particle distance dependence and an interaction coefficient ε; and(b) performing, via the processor. a molecular dynamics calculation on arenormalized particle system S′ in which the number of renormalizedparticles is N′, the mass of each particle is the mass m′ of eachrenormalized particle, and inter-particle interaction potential energyis expressed by the product εf of the dimensionless function f and theinteraction coefficient ε.
 7. The simulation method according to claim6, further comprising the step of: (c) when the position vector of acertain particle in the renormalized particle system S′ obtained bymolecular dynamics calculation in the step (b) is q′, the momentumvector of a certain particle in the renormalized particle system S′obtained by molecular dynamics calculation in the step (b) is p′, andthe velocity vector of a certain particle in the renormalized particlesystem S′ obtained by molecular dynamics calculation in the step (b) isv′, carrying out at least one of calculation for obtaining the positionvector q of a certain particle in the particle system S by a conversionequation q=q′α, calculation for obtaining the momentum vector p of acertain particle in the particle system S by a conversion equationp=p′/α^(δ/2), and calculation for obtaining the velocity vector v of acertain particle in the particle system S by a conversion equationv=v′α^(δ/2) using at least one of the first renormalization factor a andthe third renormalization factor δ determined in the step (a).
 8. Asimulation method comprising the steps of: (d) obtaining, via aprocessor, a number N′ of renormalized particles by a conversionequation N′=N/(α_(x)α_(y)α_(z)), and further determining arenormalization factor δ equal to or greater than 0 to obtain a massm_(x)′ of each renormalized particle regarding an X direction by aconversion equation m_(x)′=mα_(x) ^(δ), to obtain a mass m_(Y)′ of eachrenormalized particle regarding a Y direction by a conversion equationm_(Y)′=mα_(Y) ^(δ), and to obtain a mass m_(Z)′ of each renormalizedparticle regarding a Z direction by a conversion equation m_(Z)′=mα_(Z)⁶⁷ , wherein a particle system S in which a number of particles is N, amass of each particle is m, and inter-particle interaction potentialenergy can be expressed by a product εf of a dimensionless function frepresenting inter-particle distance dependence and an interactioncoefficient ε, and an XYZ orthogonal coordinate system in a space wherethe particle system S is arranged, on the basis of a renormalizationfactor α_(x) in an X direction, a renormalization factor α_(Y), in a Ydirection, and a renormalization factor α_(Z) in a Z direction greaterthan 1; and (e) performing, via the processor, a molecular dynamicscalculation on a renormalized particle system S′, in which the number ofrenormalized particles is N′, the mass of each particle is a mass m′ ofeach renormalized particle, and inter-particle interaction potentialenergy is expressed by the product εf of the dimensionless function fand an interaction coefficient ε, using the mass m_(x)′ of eachrenormalized particle regarding the X direction for an equation ofmotion in the X direction, the mass m_(Y)′ of each renormalized particleregarding the Y direction for an equation of motion in the Y direction,and the mass m_(Z)′ of each renormalized particle regarding the Zdirection for an equation of motion in the Z direction.
 9. Thesimulation method according to claim 8, further comprising the step of:(f) when the position vector of a certain particle in the renormalizedparticle system S′ obtained by molecular dynamics calculation in thestep (e) is q′, the components in the X, Y, and Z directions of theposition vector q′ are respectively q_(X)′, q_(Y)′, and q_(Z)′, themomentum vector of a certain particle in renormalized the particlesystem S′ obtained by molecular dynamics calculation in the step (e) isp′, the components in the X, Y, and Z directions of the momentum vectorp′ are respectively p_(X)′, p_(Y)′, and p_(Z)′, the velocity vector of acertain particle in the renormalized particle system S′ obtained bymolecular dynamics calculation in the step (e) is v′, and the componentsin the X, Y, and Z directions of the velocity vector v′ are respectivelyv_(X)′, v_(Y)′, and v_(Z)′, carrying out at least one of calculation forrespectively obtaining the components q_(X), q_(Y), and q_(Z) in the X,Y, and Z directions of the position vector q of a certain particle inthe particle system S by conversion equations q_(X)=q_(X)′α_(X),q_(Y)=q_(Y)′α_(Y), and q_(Z)=q_(Z)′α_(Z), calculation for respectivelyobtaining the components p_(X), p_(Y), and p_(Z) in the X, Y, and Zdirections of the momentum vector p of a certain particle in theparticle system S by conversion equations p_(X)=p_(X)′/α_(X) ^(δ/2),p_(Y)=p_(Y)′/α_(Y) ^(δ/2), and p_(Z)=p_(Z)′/α_(Z) ^(δ/2), andcalculation for respectively obtaining the components v_(X), v_(Y), andv_(Z) in the X, Y, and Z directions of the velocity vector v of acertain particle in the particle system S by conversion equationsv_(X)=v_(X)′α_(X) ^(δ/2), v_(Y)=v_(Y)′α_(Y) ^(δ/2), andv_(Z)=v_(Z)′α_(Z) ^(δ/2) using the renormalization factor α_(X) in the Xdirection, the renormalization factor α_(Y) in the Y direction, therenormalization factor α_(Z) in the Z direction, and the renormalizationfactor ε determined in the step (d).
 10. The simulation method accordingto claim 8, wherein at least two of the renormalization factor α_(x) inthe X direction, the renormalization factor α_(Y) in the Y direction,and the renormalization factor α_(Z) in the Z direction are differentfrom each other.